require(data.table)
library("dplyr")
library(ggplot2)
##############################
############################  Figure 3e boxplot
##############################

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"


#******************************************************************************
#=============================Statistics and Figures===========================
#******************************************************************************

#================load data to skip above steps=================
#load dif.8crop.wm prepared by Data.Fig3-4.R
load(paste0(wd1,"/data/RData/Fig3-4.globalmean-yieldchanges.Rdata")) #including all grids with crop failures to calc modeled total change in yield under different scenarios (area>1000ha)

#load MLR results with confidence intervals from Bootstrap resampling (by Data3.MLR-log-bootstrap-resampling.R)
load(paste0(wd1,"/data/RData/MLR-pergrid-log-wtm-effect-Bootstrap-1000.Rdata"))


#=========aggregate per-crop effects from resampling to global total effects===========
eff.boots.pct <-  eff.log.boots[dif.8crop.wm,
                                .( Level, T.eff=(exp(T.log)-1)*100, R.eff=(exp(R.log)-1)*100, M.eff=(exp(M.log)-1)*100,
                                   RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, P.eff=(exp(P.log)-1)*100, RH.eff=(exp(RH.log)-1)*100,
                                   co2.eff=(exp(i.co2.log)-1)*100, luc.eff=(exp(i.luc.log)-1)*100, #note: i.luc.log from dif.8crop (all based on SSP5 surface data) is still not the total luc effect because SSP2 have different crop grid cells than SSP5
                                   Tot.eff=(exp(Tot.log)-1)*100, Y.mod.dif=(exp(i.Y.mod.dif)-1)*100, 
                                   Tot.eff2=(exp(Tot.log2)-1)*100, Y.mod.dif2=(exp(i.Y.mod.dif2)-1)*100, #wtm of per-grid difference i.Y.mod.dif is higher than differene of global wtm yield (yield-yield0)/yield0*100
                                   Tot.eff3=(exp(Tot.log3)-1)*100, #per grid luc.log from dif.8crop.s included in tot.log3 (8.4%), larger than adding per crop i.luc.log from dif.8crop to tot.eff3 (7.4%)
                                   Y.mod.dif3=(exp(i.Y.mod.dif3)-1)*100, #per grid wtm modeled change from dif.8crop, slightly higher than (exp(i.Y.mod.dif+i.co2.log+i.luc.log)-1)*100
                                   yield=i.yield, yield1=i.yield1, yield0=i.yield0, area=i.area ,yield2=i.yield2, area2=i.area2 #prod=i.prod, prod0=i.prod0, prod1=i.prod1,
                                ), by=.EACHI, #on=.(Scenario,Year)]
                                on=.(Scenario,Crop,Irrig,Year)]
eff.boots.pct <- eff.boots.pct[!is.na(Level)]
#aggregate per crop effect to global considering base yield difference between crop types
eff.tot <- rbind(eff.boots.pct[Scenario=="RCP45 - RCP85",
                               .(T.eff=weighted.mean(T.eff, yield1*area, na.rm = T), R.eff=weighted.mean(R.eff, yield1*area, na.rm = T), M.eff=weighted.mean(M.eff, yield1*area, na.rm = T),
                                 RD.eff=weighted.mean(RD.eff, yield1*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield1*area, na.rm = T),
                                 P.eff=weighted.mean(P.eff, yield1*area, na.rm = T), RH.eff=weighted.mean(RH.eff, yield1*area, na.rm = T),
                                 Tot.eff=weighted.mean(Tot.eff, yield1*area, na.rm = T), luc.eff=weighted.mean(luc.eff, yield*area, na.rm = T), #LUC has differnet base yield from RCP45_85LU
                                 Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), co2.eff=weighted.mean(co2.eff, yield0*area, na.rm = T),
                                 Tot.eff3=weighted.mean(Tot.eff3, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield1*area, na.rm = T),
                                 Y.mod.dif2=weighted.mean(Y.mod.dif2, yield0*area, na.rm = T), Y.mod.dif3=weighted.mean(Y.mod.dif3, yield0*area, na.rm = T),
                                 yield1=weighted.mean(yield1, area), yield0=weighted.mean(yield0, area), yield=weighted.mean(yield, area), 
                                 yield2=weighted.mean(yield2, area2), area2=sum(area2, na.rm=T),
                                 area=sum(area, na.rm=T)), by =.(Scenario,Level,Year)],
                 eff.boots.pct[Scenario!="RCP45 - RCP85",
                               .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),M.eff=weighted.mean(M.eff, yield0*area, na.rm = T),
                                 RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),
                                 P.eff=weighted.mean(P.eff, yield0*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield0*area, na.rm = T),
                                 Tot.eff=weighted.mean(Tot.eff, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield0*area, na.rm = T),
                                 yield0=weighted.mean(yield0, area), yield=weighted.mean(yield, area), luc.eff=0, co2.eff=0,
                                 area=sum(area, na.rm=T)),
                               by =.(Scenario,Level,Year)], fill=TRUE)



#======Version 1: only show Climate+CO2 effect for ER (RCP45(SSP5LU))
eff.tot3 <- copy(eff.tot)
eff.tot3 <- eff.tot3[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff2] #statistical estimate of grid-level wtm effect
eff.tot3 <- eff.tot3[, Y.mod.dif:=(yield-yield0)/yield0*100] #relative change of global average yield (yield=RCP45_85LU,SAI,MSB,CCT, yield0=RCP85)

#=======statistics based on Fig.3e
eff.tot3[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff), by=.(Scenario)]  #climate only
eff.tot3[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff2), by=.(Scenario)] #climate+CO2 for ER
eff.tot3[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff3), by=.(Scenario)] #climate+CO2+LUC for ER
eff.tot3[Year >=2075 & Year<=2099 & Level=="Mean", mean(co2.eff), by=.(Scenario)]
eff.tot3[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] #CLM5 simulated total effect
#===impact on global yield: 11%, 9%, 11%, and -5% for SAI, MSB, CCT and ER, respectively

#report the error term in Figure 3 (average cross all crops)
eff.tot3[Year >=2075 & Year<=2099 & Level=="Mean", .(Residual=mean(Tot.eff-Y.mod.dif)), by=.(Scenario)] #difference of statistical estimate and modeled change
#===residual errors: -0.8%, 1.1%, 0.3%, -0.6% for SAI, MSB, CCT and ER, respectively


#% effect on global yield
eff.tot3.t <- eff.tot3[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Year,Scenario)]
names(eff.tot3.t)[3:6] <- c("Variable", "Mean", "Conf.L",  "Conf.U")

casecolor1 <- c("#CC79A7", "#8491B4FF", "#91D1C2FF","#7E6148FF")
barplot(height=rep(1,length(casecolor1)),col=casecolor1) #display colors

levs2 <- c("T.eff",  "RD.eff", "RI.eff", "P.eff", "RH.eff", "co2.eff", "Tot.eff", "Y.mod.dif") #, "P.mod.dif")#
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle(' Total '), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle(' Total '), textstyle('(CLM5)')), ' '))) # extstyle tells not to shrink the text size when atop more than 2 layers

levs3 <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85","RCP45 - RCP85")
scenario <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "Emissions reduction")


df <-  eff.tot3.t[Variable%in%levs2 & Year >=2075 & Year<=2099]
p <- df %>%  #boxplot show the percentiles of 2075-2099 mean values
  mutate(Scenario = factor(Scenario, levels=levs3, labels=levs3)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  #ggplot(aes(x=Crop, y=Mean, fill=Scenario)) +
  ggplot(aes(x=Variable, y=Mean, fill=Scenario)) +
  geom_boxplot(notch = FALSE, lwd=0.3) +
  geom_hline(yintercept=0) + geom_vline(aes(xintercept=6.5),linetype="dashed") +
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  scale_y_continuous("Change in global yield (%)")+#, limits=c(-25, 25), breaks = c(-20, -10, 0, 10, 20, 30)) +
  scale_fill_manual(name="", values= casecolor1, labels=scenario) +
  #scale_fill_manual(name="", values= rev(c("#26828EFF","#B4DE2CFF", "#56B4E9","#CC6600")), labels=scenario) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=1, keywidth=1)) #nrow =1))
p + #coord_flip() +
  theme_bw()+ #theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression(''*bold(e)*'')) + #"(b) Partial and total climate effects on rainfed crops in 2090-2099"
  theme(plot.title = element_text(size = 16,hjust=0.), legend.text = element_text(size=15), legend.position = c(0.5, -0.3)) + #use numeric vector c(x,y)
  theme(axis.title = element_text(size=15), strip.text.x = element_text(size = 15, vjust=0)) +
  theme(axis.text = element_text(size=15), axis.text.x = element_text(angle=0,hjust=0.5,vjust=0.5))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"))
#use specific legend.position = c(0.5, -0.3) to adjust the location of legend at any place
#c(0,0) corresponds to the “bottom left” and c(1,1) corresponds to the “top right” position

ggsave(paste0(wd1,"/figures/main/Fig3e.svg"), device ="svg", units = "in", width = 8, height = 4.2)




##save as PDF (1-column, 8.8cm)--resize line thickness and text sizes
df <-  eff.tot3.t[Variable%in%levs2 & Year >=2075 & Year<=2099]
p <- df %>%  #boxplot show the percentiles of 2075-2099 mean values
  mutate(Scenario = factor(Scenario, levels=levs3, labels=levs3)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  ggplot(aes(x=Variable, y=Mean, fill=Scenario)) +
  geom_boxplot(notch = FALSE, lwd=0.1, outlier.size = 0.3) +
  geom_hline(yintercept=0, size=0.3) + geom_vline(aes(xintercept=6.5),linetype="dashed", size=0.3) +
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  scale_y_continuous("Change in global yield (%)")+#, limits=c(-25, 25), breaks = c(-20, -10, 0, 10, 20, 30)) +
  scale_fill_manual(name="", values= casecolor1, labels=scenario) +
  #scale_fill_manual(name="", values= rev(c("#26828EFF","#B4DE2CFF", "#56B4E9","#CC6600")), labels=scenario) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=0.5, keywidth=0.5)) #nrow =1))

p + #coord_flip() +
  theme_bw(base_size = 5)+ theme(element_line(size = 0.2)) + #theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression(''*bold(e)*'')) + #"(b) Partial and total climate effects on rainfed crops in 2090-2099"
  theme(plot.title = element_text(size = 7,hjust=0), 
        legend.text = element_text(size=6), 
        legend.position = c(0.5, -0.3)) + #use numeric vector c(x,y)
  theme(axis.title = element_text(size=7)) +
  theme(axis.text = element_text(size=7), axis.text.x = element_text(angle=0,hjust=0.5,vjust=0.5))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"))

ggsave(paste0(wd1,"/figures/main/Figure3e.svg"), device ="svg", units = "in", width = 3.5, height = 2)
ggsave(paste0(wd1,"/figures/main/Figure3e.pdf"), device ="pdf", units = "in", width = 3.5, height = 2)






#========================normalizing/scaling by TOA radiative imbalalnce===================
TOM <- read.table(paste0(wd1,"/data/NCL/TOM_radiative_imbalance.2020-2099.txt"), header = TRUE, sep=",")
TOM <- as.data.table(TOM)
TOM$Year <- 2020:2099
TOM.f <- TOM[,.(ER=RCP45-RCP85, SAI=RCP85_SAI-RCP85, MSB=RCP85_MSB-RCP85, CCT=RCP85_CCT-RCP85), by=.(Year)]
#ER,SAI,MSB all have negative values, only CCT has a few positive values (RTMT > RCP85)
#original TOM radiative flux imbalance (not the same as sustained radiative forcing in equilibrum)
plot(ER~Year, data=TOM.f, type="l", col= "#26828EFF", ylim=c(-2,0.2))
lines(SAI~Year, TOM.f, col="#CC6600")
lines(MSB~Year, TOM.f, col="#56B4E9")
lines(CCT~Year, TOM.f, col="#B4DE2CFF")

#==Method: rescale the net effects by matching the mean and trend of SAI/MSB/CCT to those of ER
TOM.t <- TOM.f[, data.table(t(.SD), keep.rownames=TRUE), by=.(Year)]
names(TOM.t)[2:3] <- c("Scenario", "Data")

TOM.lm <- TOM.t[, .(intercept=lm(Data~Year)$coefficients[1], slope=lm(Data~Year)$coefficients[2],
                    Mean=mean(Data, na.rm = T)), by =.(Scenario)]
TOM.dtr <- TOM.t[TOM.lm, .(Data.dtr1=Data-i.slope*Year-i.intercept+i.Mean, #only remove the trend but keep the original mean
                           Data.dtr2=Data-i.slope*Year-i.intercept, #same as detrend() function, mean becomes zero
                           Data, intercept=i.intercept,slope=i.slope,Mean=i.Mean, Year), by=.EACHI, on =.(Scenario)]

#adjust the mean and the rescale the slope of SAI/MSB/CCT to those of ER
ER.lm <- TOM.lm[Scenario=="ER",] #baseline
TOM.scaler <- TOM.lm[, .(mean.adj = ER.lm$Mean-Mean, slope.adj = ER.lm$slope-slope, 
                         mean.ratio = ER.lm$Mean/Mean, slope.ratio = ER.lm$slope/slope), by=.(Scenario)]
#scaling method 1: return detrended data to scaled values
TOM.scaled1 <- TOM.dtr[TOM.scaler, Data.new1:=Data.dtr2+slope*Year*slope.ratio + ER.lm$intercept, #add the intercept of ER
                      by=.EACHI, on=.(Scenario)]
#scaling method 2: avoid using intercept of lm, only need Year and original Data, mean(Year)=2059.5 (2020:2099)
TOM.scaled2 <- TOM.dtr[TOM.scaler, Data.new2:=Data+mean.adj+(Year-mean(Year))*slope.adj, #same as above method
                      by=.EACHI, on=.(Scenario)]
#but mean.adj is specific to the TOM data, cannot be applied to the new data eff.tot
#scaling method 3: only need Year and original Data
TOM.scaled3 <- TOM.dtr[TOM.scaler, Data.new3:=Data+mean(Data,na.rm = T)*(mean.ratio-1)+(Year-mean(Year))*slope.adj, #same as above method
                      by=.EACHI, on=.(Scenario)]
#method 3 is the best way to use (do not need to do lm regression and can directly rescale the mean and slope)
#however, slope.adj is still specific to the TOM data, cannot be applied to the new data eff.tot

#now all scenarios have the same mean and slope 
TOM.scaled1[, .(Mean=mean(Data.new1, na.rm=T), slope=lm(Data.new1~Year)$coefficients[2]), by=.(Scenario)]
TOM.scaled2[, .(Mean=mean(Data.new2, na.rm=T), slope=lm(Data.new2~Year)$coefficients[2]), by=.(Scenario)]
TOM.scaled3[, .(Mean=mean(Data.new3, na.rm=T), slope=lm(Data.new3~Year)$coefficients[2]), by=.(Scenario)]
#(mean.ratio-1) can be applied to the mean of new data (eff.tot), 
#but slope.adj cannot be directly applied to partial and total effects because the slope of each effect variable is different from that of TOA imbalance

#plot rescaled time series of TOM radiative imbalance change relative to RCP85
plot(ER~Year, data=TOM.f, type="l", col= "#26828EFF", cex.main=1.0,
     main="TOA imbalance \n difference from RCP85", ylim=c(-2,0.2), ylab=expression("W m"^-2))
lines(SAI~Year, TOM.f, col="#CC6600", lwd=1.5)
lines(MSB~Year, TOM.f, col="#56B4E9", lwd=1.5)
lines(CCT~Year, TOM.f, col="#B4DE2CFF", lwd=1.5)
lines(Data.new2~Year, TOM.scaled3[Scenario=="SAI"], col="#CC6600", lty=2, lwd=2)  
lines(Data.new2~Year, TOM.scaled3[Scenario=="MSB"], col="#56B4E9", lty=3, lwd=2)   
lines(Data.new2~Year, TOM.scaled3[Scenario=="CCT"], col="#B4DE2CFF", lty=4, lwd=2)  


#======Finally use the ratio between the fit line of SAI/MSB/CCT and the fit line of ER as single scaler on any effect variable
eff.scaler <- TOM.scaled3[, .(Year, scaler= lm(Data.new3~Year)$fitted.values/lm(Data~Year)$fitted.values,
                           fit=lm(Data~Year)$fitted.values, fit.new=lm(Data.new3~Year)$fitted.values), by=.(Scenario)] 
#the "scaler" scales the RTMT(TOA imbalance) mean and slope (i.e. the fit line) of SGs to that of ER
#all scenarios now have the identical fit line (same fit.new)
lines(fit.new~Year, eff.scaler[Scenario=="ER"], col="#26828EFF", lty=1, lwd=2)
lines(fit~Year, eff.scaler[Scenario=="SAI"], col="#CC6600", lty=1, lwd=2)  
lines(fit~Year, eff.scaler[Scenario=="MSB"], col="#56B4E9", lty=1, lwd=2)   
lines(fit~Year, eff.scaler[Scenario=="CCT"], col="#B4DE2CFF", lty=1, lwd=2)  
lines(fit.new~Year, eff.scaler[Scenario=="SAI"], col="#CC6600", lty=2, lwd=2)  
lines(fit.new~Year, eff.scaler[Scenario=="MSB"], col="#56B4E9", lty=3, lwd=2)   
lines(fit.new~Year, eff.scaler[Scenario=="CCT"], col="#B4DE2CFF", lty=4, lwd=2)  


#===apply scaler to eff.tot3  
eff.scaler[Scenario=="ER",Scenario:="RCP45 - RCP85"]
eff.scaler[Scenario=="SAI",Scenario:="SAI - RCP85"]
eff.scaler[Scenario=="MSB",Scenario:="MSB - RCP85"]
eff.scaler[Scenario=="CCT",Scenario:="CCT - RCP85"]

eff.tot3.scaled <- copy(eff.tot3.t)
eff.tot3.scaled <- eff.tot3.scaled[eff.scaler, scaler:=i.scaler, by=.EACHI, on=.(Scenario)]

#==apply the scaler to all partial and total effects including CLM modeled difference, but excluding base yield and crop area
eff.tot3.scaled <- eff.tot3.scaled[!Variable%in%c("yield1", "yield","yield2","area2","yield0","area"),
                                   .(Year, Mean=Mean*scaler, Conf.L=Conf.L*scaler, Conf.U=Conf.U*scaler),
                                   by=.(Scenario,Variable)]

#double check the scaling factor is consistent at any year for a given Scenario
tail(eff.scaler[Scenario=="CCT - RCP85"])
TOM.scaled3[Year==2095 & Scenario=="CCT",Data.new3]/TOM.scaled3[Year==2095 & Scenario=="CCT",Data] #1.42
eff.scaler[Year==2095 & Scenario=="CCT - RCP85",fit.new]/eff.scaler[Year==2095 & Scenario=="CCT - RCP85",fit] #1.46
eff.tot3.scaled[Year==2095 & Scenario=="CCT - RCP85" & Variable=="Tot.eff",Mean]/eff.tot3.t[Year==2095 & Scenario=="CCT - RCP85" & Variable=="Tot.eff", Mean] #1.46
eff.tot3.scaled[Year==2095 & Scenario=="CCT - RCP85" & Variable=="T.eff", Mean]/eff.tot3.t[Year==2095 & Scenario=="CCT - RCP85" & Variable=="T.eff", Mean] #1.46

eff.tot3.new <- rbind(eff.tot3.scaled, eff.tot3.t[Variable%in%c("yield1", "yield","yield2","area2","yield0","area")])



levs2 <- c("T.eff",  "RD.eff", "RI.eff", "P.eff", "RH.eff", "co2.eff", "Tot.eff", "Y.mod.dif") #, "P.mod.dif")#
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle(' Total '), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle(' Total '), textstyle('(CLM5)')), ' '))) # extstyle tells not to shrink the text size when atop more than 2 layers

levs3 <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85","RCP45 - RCP85")
scenario <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "Emissions reduction")

df <-  eff.tot3.new[Variable%in%levs2 & Year >=2075 & Year<=2099]
p <- df %>%  #boxplot show the percentiles of 2075-2099 mean values
  mutate(Scenario = factor(Scenario, levels=levs3, labels=levs3)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  #ggplot(aes(x=Crop, y=Mean, fill=Scenario)) +
  ggplot(aes(x=Variable, y=Mean, fill=Scenario)) +
  geom_boxplot(notch = FALSE, lwd=0.3) +
  geom_hline(yintercept=0) + geom_vline(aes(xintercept=6.5),linetype="dashed") +
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  scale_y_continuous("Change in global yield (%)")+#, limits=c(-25, 25), breaks = c(-20, -10, 0, 10, 20, 30)) +
  scale_fill_manual(name="", values= casecolor1, labels=scenario) +
  #scale_fill_manual(name="", values= rev(c("#26828EFF","#B4DE2CFF", "#56B4E9","#CC6600")), labels=scenario) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=1, keywidth=1)) #nrow =1))
p + #coord_flip() +
  theme_bw()+ #theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression('('*bold(e)*')')) + #"(b) Partial and total climate effects on rainfed crops in 2090-2099"
  theme(plot.title = element_text(size = 16,hjust=0.5), legend.text = element_text(size=15), legend.position = c(0.5, -0.3)) + #use numeric vector c(x,y)
  theme(axis.title = element_text(size=16), strip.text.x = element_text(size = 16, vjust=0)) +
  theme(axis.text = element_text(size=16), axis.text.x = element_text(angle=0,hjust=0.5,vjust=0.5))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_rect(colour = "black"))

#Fig.3e scaled by TOA radiative imbalance
ggsave(paste0(wd1,"/figures/supplementary/Fig.S11e.svg"), device ="svg", units = "in", width = 8, height = 4.2)




#===(Supplement Fig. S3a) all effects climate+CO2+LUC for ER
eff.tot2 <- copy(eff.tot)
eff.tot2 <- eff.tot2[Scenario=="RCP45 - RCP85", Tot.eff:=Tot.eff3] #statistical estimate of grid-level wtm effect
eff.tot2 <- eff.tot2[Scenario=="RCP45 - RCP85", Y.mod.dif:=(yield2-yield0)/yield0*100] #including CO2 and LUC effects for ER
eff.tot2 <- eff.tot2[Scenario!="RCP45 - RCP85", Y.mod.dif:=(yield-yield0)/yield0*100]

#add total effect on global production (considering all crop grid cells of SSP2 and SSP5)
eff.tot2 <- eff.tot2[Scenario=="RCP45 - RCP85", P.mod.dif:=(yield2*area2-yield0*area)/(yield0*area)*100]
eff.tot2 <- eff.tot2[Scenario!="RCP45 - RCP85", P.mod.dif:=(yield*area-yield0*area)/(yield0*area)*100]
#add residual effects including residual errors in the regression and changes of samples not included (those exist <60 years)
eff.tot2 <- eff.tot2[, Residual:=Y.mod.dif-Tot.eff] #difference of statistical estimate and modeled change

#% effect on global yield
eff.tot2.t <- eff.tot2[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Year,Scenario)]
names(eff.tot2.t)[3:6] <- c("Variable", "Mean", "Conf.L",  "Conf.U")
#names(eff.tot2.t)[3:4] <- c("Variable", "Mean")

eff.tot2[Year >=2075 & Year<=2099 & Level=="Mean", mean(Tot.eff), by=.(Scenario)]
eff.tot2[Year >=2075 & Year<=2099 & Level=="Mean", mean(luc.eff), by=.(Scenario)]
eff.tot2[Year >=2075 & Year<=2099 & Level=="Mean", mean(Y.mod.dif), by=.(Scenario)] #6-11% on yield
eff.tot2[Year >=2075 & Year<=2099 & Level=="Mean", mean(P.mod.dif), by=.(Scenario)] #6-11% on production


#======Boxplot for global effects (climate+co2+LUC)
levs2 <- c("T.eff",  "RD.eff", "RI.eff", "P.eff", "RH.eff", "co2.eff", "luc.eff", "Tot.eff", "Y.mod.dif")
varname <- c(expression(atop(atop(textstyle('T'), textstyle('effect')), ' ')), #expression(atop("Diffuse \n radiation", 'effect')),
             expression(atop(atop(textstyle('RD'), textstyle('effect')),' ')),
             expression(atop(atop(textstyle('RI'),textstyle('effect')),' ')),
             expression(atop(atop(textstyle('P'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('RH'), textstyle('effect')),  ' ')),
             expression(atop(atop(textstyle('CO'[2]), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('LUC'), textstyle('effect')), ' ')),
             expression(atop(atop(textstyle('Total effect'), textstyle('(MLR)')), ' ')),
             expression(atop(atop(textstyle('Total effect'), textstyle('(CLM5)')), ' '))) # textstyle tells not to shrink the text size when atop more than 2 layers
levs3 <- c("SAI - RCP85", "MSB - RCP85", "CCT - RCP85", "RCP45 - RCP85")
#scenario <- c( "SAI - RCP85", "MSB - RCP85", "CCT - RCP85" ,"Emissions reduction")
scenario = levs3

df <-  eff.tot2.t[Variable%in%levs2 & Year >=2075 & Year<=2099]
p <- df %>%  #boxplot show the percentiles of 2075-2099 mean values
  mutate(Scenario = factor(Scenario, levels=levs3, labels=scenario)) %>%
  mutate(Variable = factor(Variable, levels=levs2, labels=varname )) %>%
  #ggplot(aes(x=Crop, y=Mean, fill=Scenario)) +
  ggplot(aes(x=Variable, y=Mean, fill=Scenario)) +
  geom_boxplot(notch = FALSE) +
  geom_hline(yintercept=0) + geom_vline(aes(xintercept=7.5),linetype="dashed") +
  scale_x_discrete(name="", labels=varname) + #label_parsed only works with facet_wrap or facet_grid
  scale_y_continuous("Change in global yield (%)", limits=c(-25, 35), breaks = c(-20, -10, 0, 10, 20, 30)) +
  #scale_fill_manual(name="", values= rev(c("#26828EFF","#B4DE2CFF", "#56B4E9","#CC6600")), labels=scenario) +
  scale_fill_manual(name="", values= casecolor1, labels=scenario) +
  guides(fill = guide_legend(title=NULL,nrow=1,byrow=T,label.position = "right", label.hjust = 0, keyheight=1, keywidth=1)) #nrow =1))
p + #coord_flip() +
  theme_bw()+ #theme_minimal()+  #theme_classic() +  #horizontal bars
  ggtitle(expression('('*bold(a)*')')) + #"(b) Partial and total climate effects on rainfed crops in 2090-2099"
  theme(plot.title = element_text(size = 18,hjust=0.5)) +
  theme(axis.title = element_text(size=18), strip.text.x = element_text(size = 18, vjust=0)) +
  theme(axis.text = element_text(size=18), axis.text.x = element_text(angle=0,hjust=0.5,vjust=1))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.background = element_blank(), #legend.position = "bottom"
        legend.position = c(0.35,0.9), legend.text = element_text(size=18),
        panel.border = element_rect(colour = "black"))

ggsave(paste0(wd1,"/figures/supplementary/Fig.S3a.svg"), device ="svg", units = "in", width = 13, height = 5.5)






